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We describe a systematic approach for the efficient numerical solution of nonlinear Schrodinger- 
type partial differential equations of the form (K, + V + g\ip\ 2 )ip = 0, with an energy operator K., 
a scalar potential V, and a scalar parameter g. Instrumental to the approach are developments in 
numerical linear and nonlinear algebra, specifically numerical parameter continuation. We demon- 
strate how a continuous sequence of solutions can be obtained regardless of their stability, so that 
finally the spectrum of stable and unstable solutions in the specified parameter range is fully re- 
vealed. The method is demonstrated for the GL equation in a three-dimensional superconducting 
domain with an inhomogeneous magnetic field, a numerically demanding problem known to have 
an involved solution landscape. 
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Nonlinear Schrodinger equations and their variations 
are used to model a wide variety of physical systems ^lj.2], 
with applications spanning superconductivity [3 , quan- 
tum condensates [H [5] , nonlinear acoustics [5] , nonlinear 
optics [7| , and hydrodynamics 8 . Typically, the physical 
models described by these equations contain a set of pa- 
rameters specifying, e.g., the sample geometry, external 
fields, or boundary conditions. To understand the phys- 
ical properties of the system, it is interesting to explore 
the energy landscape of the steady state solutions as a 
function of one or more of these parameters. In general, 
slight perturbations in one of the control parameters can 
induce changes in the stability properties of the states, 
causing abrupt transitions in the energy landscape [3]. 
Although the existence of steady states can be proven in 
certain cases, analytic solutions are hard or impossible to 
obtain in realistic systems. Numerical methods are hence 
of particular importance for understanding the physics of 
the systems modeled by nonlinear Schrodinger equations. 
The challenge is to develop efficient computational tools 
to explore the full energy landscape, including minima 
and saddle-points, of three- and higher-dimensional sys- 
tems. 

In general, the nonlinear Schrodinger equations de- 
scribing the evolution of a quantum-dynamical system 
represented by a complex-valued order parameter vb are 
written as 

with a linear, positive-semidefinite (energy) operator K., 
an external potential V, and the coupling parameter g. 
The term \ip\ 2 usually describes a probability density of 
the model entity, e.g., the locality of quantum particles. 
Examples include the Gross-Pitaevskii equation where 
K, = — A is the negative Laplacian, and the Ginzburg- 



Landau (GL) equation, where K, = (— iV — A) 2 is the 
covariant Laplacian with a given vector potential A. 

To understand the long-term dynamics, it is essential 
to compute steady states of the system, i.e., solutions to 

= S(4>). (2) 

Solving ([2]) in realistic three-dimensional domains is a dif- 
ficult numerical task: the number of unknowns quickly 
becomes very large and standard numerical methods be- 
come impractical. In addition, the energy landscape 
of the solutions becomes very complicated and its sys- 
tematic exploration is prohibitive with current numerical 
techniques. 

This letter describes a numerical approach for the ef- 
ficient computation of the steady-state landscape as a 
function of the control parameters. Its central com- 
ponent is a Newton-Krylov algorithm [TDJ HI] to solve 
the nonlinear problem (pi), combined with numerical pa- 
rameter continuation [12] to explore the solution land- 
scape. Although numerical continuation and Newton- 
Krylov solvers are well-known methods for large-scale 
systems, they cannot be applied straightforwar dly to ((2|. 
We will show that, by exploiting the properties of the lin- 
earization of the operator S(-) and devising a specially 
tailored preconditioner, it is possible to considerably ac- 
celerate the convergence of the linear iterations. This 
opens up the possibility to compute steady states and 
to systematically explore the energy landscape in three- 
dimensional problems. The new method scales optimally 
with the number of unknowns in the system and is fully 
parallclizable. In particular, we illustrate the power of 
the method by studying three-dimensional vortex nucle- 
ation in an extreme-type-II superconductor in an inho- 
mogeneous magnetic field, described by the GL equation. 

Numerical simulations within the GL model are an es- 
sential tool for the analysis of superconducting phenom- 
ena. In this area, vortex matter has been at the forefront 
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of research in the past two decades. Emphasis has been 
put on the computation of vortex states with imposed 
confinement (i.e., the sample shape) and on their depen- 
dence upon critical parameters of the superconducting 
sample. Of particular relevance are unstable states of 
the system, often called saddle points. These solutions 
shape the energy landscape as they constitute the connec- 
tions between families of stable states, thereby providing 
a unique insight into the dynamic transitions and vor- 
tex rearrangements that have been observed experimen- 
tally. Notably, with the framework proposed in this letter 
we can compute both stable and unstable states, which 
are not accessible with traditional numerical methods. 
Owing to numerical difficulties, saddle points have been 
calculated only for radially-symmetric samples such as 
disks, using a limited-expansion method |13j . Recently, 
the full energy landscape (including saddle points) was 
systematically explored in two-dimensional square sam- 
ples, using numerical continuation techniques [14j : in 
particular, it was possible to build an atlas of the instabil- 
ities occurring in the sample, providing a complete clas- 
sification of the symmetries of observable stable states. 
In this letter, we address the much more challenging case 
of three-dimensional samples of arbitrary shape. 

Existing methods The literature on numerical meth- 
ods for the time-dependent equation (JlJ is rather exten- 
sive and mostly concerned with time-stepping schemes 
[TMT7] . For example, references [THJ [TU] leveraged 
specific properties of certain numerical procedures and 
settings for dealing with the Gross-Pitaevskii equa- 
tion. Stationary-states are typically found by applying 
a (pseudo-)time-stepping scheme until a stationary state 
is reached [20-22 . There are, however, several disad- 
vantages with this approach. Firstly, iterations converge 
only for strictly stable states, therefore unstable or saddle 
point states can not be computed. Secondly, stable solu- 
tions may have extremely long (sometimes oscillatory) 
transients, therefore convergence for three-dimensional 
domains may be prohibitively slow. 

Newton's method A better approach is to use New- 
ton iterations directly on |2]), starting from a suitable 
initial guess: Newton's method converges superlinearly 
in a neighborhood of the solution, irrespectively of the 
stability properties of the equilibrium. Once a steady 
state ip is found, stability is determined by computing 
the spectrum of the operator obtained by linearizing S 
around For the Schrodinger equations, this linear op- 
erator is defined via the action 

J{^={K, + V + 2g\^\ 2 )^ + g^ (3) 

where <f> denotes complex conjugation. 

A sequence of approximations ip^ to the steady state 
is computed with Newton's method = ip^ + 

5ip' k \ where the update Sip^ satisfies 

J(if>W)6if>W = -S(^ k) ). (4) 



Therefore, the solution of a large linear system is required 
at each Newton step k, which is the most significant 
difficulty when applying Newton's method to nonlinear 
Schrodinger equations. 

Solving the Jacobian system Linear systems such 
as Q can be solved using Krylov iterative methods and 
have been widely used in the past decades [531 HI] ■ A 
property of Krylov subspace methods is that no explicit 
(matrix) representation of the operator is needed, but 
only its application to vectors (cf. p])). The convergence 
of those methods is highly dependent on the spectrum 
of the involved linear operator. Principal optimizations 
can be employed if all eigenvalues of the respective lin- 
ear operator are real-valued. This is the case if the lin- 
ear operator is represented by a Hermitian matrix or, in 
general, if J is self-adjoint with respect to a given inner 
product. The linear operator ^ associated with the non- 
linear Schrodinger equations is self-adjoint with respect 
to the inner product 

(^):=&Uj^J. (5) 

This suggests the use of MINRES [24i, a Krylov sub- 
space method suitable for indefinite self-adjoint prob- 
lems. However, one characteristic of Krylov methods is 
that a larger number of unknowns increases the num- 
ber of iterations that are needed to achieve convergence. 
In addition, the computational cost of a single evalu- 
ation of the linear action also grows with the number 
of unknowns. Therefore, high-resolution discretizations 
of three-dimensional systems would require a prohibitive 
computational effort. Indeed, decreasing the number of 
Krylov iterations is the subject of extensive research ef- 
forts in this area. 

A popular approach is to use a preconditioner for 
the linear problem. The main idea is that, instead of 
solving the discretized version Jx = b of @, one can 
solve an equivalent, numerically more favorable problem 
P~ x Jx = P~ l b with a linear, invertible preconditioning 
operator P. If P is appropriately chosen, Krylov methods 
applied to the new system converge much faster. In the 
case of the linearization of nonlinear Schrodinger opera- 
tors ([3]), the energy operator K is of particular interest, 
as it typically dominates the spectral behavior of J[^). 
More precisely, we define the symmetric preconditioning 
operator 

7>(V) :=/C + 2max{ 5 , £ }|V| 2 , (6) 

with < e <g 1 [25]. We note that V(ip) is positive- 
definite except for the uninteresting case of ip = 0. This, 
most importantly, makes the inversion of the discretized 
V(ip), P _1 (-0), computationally cheap since its positive- 
definiteness makes it a suitable target for geometric or al- 
gebraic multigrid (AMG) solvers that yield optimal con- 
vergence |26j . As will be shown, even inexact inversions 
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of ([6]) used as preconditioners for (|3]) make the Krylov 
convergence independent of the number of unknowns. 

Numerical parameter continuation The efficient lin- 
ear solver outlined above is an essential building block 
for the exploration of the energy landscape, which is 
performed via numerical parameter continuation, a well- 
established technique for numerical bifurcation analysis 
of dynamical systems [27]. Let F(ip;p) — be a nonlin- 
ear system dependent upon a scalar parameter p G K 
and let be a solution for a given parameter value 
Pq. Under generic regularity conditions for F, it is pos- 
sible to construct a one-parameter family of solutions ip, 
parametrized by p, in the neighborhood of {ipo,Po) [12] - 
First, a prediction step is taken in the tangent 

direction to the one-parameter family, then a correction 
is done using Newton's method. This leads to a new so- 
lution (t/>i,Pi). The set of points (tpl,pl) form a smooth 
solution curve. Once again, we remark that the method 
is oblivious to the stability properties of the solutions. 

Application to the Ginzburg-Landau problem We il- 
lustrate the power of this method on a numerically chal- 
lenging problem: the computation of vortex structures 
in a three-dimensional superconducting domain with a 
magnetic core, which establishes an internal and inho- 
mogeneous magnetic held. 

Given a bounded superconducting domain f2, the GL 
equations 

= (-iV-A) 2 V-^(l-|V| 2 ) intt, (7) 
= n • (— iV — A)ip on the surface d£l. 

describe stationary states of an extreme-type-II super- 
conductor subject to a magnetic held associated with the 
vector potential A [20] HH] • The equations are presented 
in dimensionless form: distances are scaled by the su- 
perconducting coherence length £, the order parameter 
-0 by its value in the absence of applied magnetic held, 
and the vector potential by £,H C 2 , where H C 2 denotes the 
upper critical magnetic held of a bulk material. Since the 
kinetic energy operator (— iV— A) 2 is Hermitian and pos- 
itive semi-definite, equation Q is of the form ^ (with 
V = — 1, g = 1) and can be solved with our numerical 
method. 

We choose f2 to be a cube with side length 10 and a 
spherical cavity of radius 1 containing a magnetic dipole 
with magnetic moment m = (0, 0, 1) T (see figure]!]). The 
associated magnetic vector potential of such a dipole is 
Aa{x) := ||a;|| _3 (m x x). This example is of great rel- 
evance to the field, since the expected loops appear in 
several physical systems [25] • Their nucleation, growth, 
motion, and recombination harbors a vast variety of 
novel physics. We perform numerical experiments us- 
ing a finite- volume (tetrahedral) discretization where the 
complex-valued order parameter -0 is approximated in 
the grid nodes [25] . 

The first important result is shown in figure [2J and 
concerns the efficiency of solving the Jacobian systems 




FIG. 1: Cube with a spherical cavity, representing a 
superconducting sample hosting a magnetic dipole 
(clipped display). The field lines of the magnetic field 
V x Ad of the dipole are shown in white. 
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Number of unknowns 

FIG. 2: The number of MINRES iterations necessary 
for converging the initial guess tp = for the system 
J(il))ip = b, ip = 1, b = 1 to a relative residual of 10~ 10 
for different preconditioners. The figure shows that the 
number of iterations increases without a preconditioner 
| o | , and remains constant for preconditioners with 
exact d h as well as approximate inversion I » |. 



with preconditioners based on the discretization P{ip) of 
More specifically, we use two preconditioners, the 
first one being the exact inverse (up to machine preci- 
sion) of P(ip) and the second one being an approximate 
inverse of P(ip) obtained with just a single AMG step. 
Preconditioned linear systems are solved for increasing 
number of unknowns and they are compared to the case 
without preconditioners. A remarkable result is that, in 
both preconditioned cases, the number of iterations does 
not increase with the number of unknowns in the system. 
This indicates optimal scalability of the solver, which is 
extremely advantageous compared to the case without 
preconditioner. 

Numerical parameter continuation is then applied to 
a discretization with 0.4 x 10 6 grid points. The mag- 



4 



netic vector potential Ad (and thus the corresponding 
magnetic field V x Ad) is scaled with the dimension- 
less magnetic moment /i = rn/(H C 2^ 3 ) which is taken as 
control parameter. For [i = 0, the homogeneous state 
ip = 1 is clearly a solution (independently of the domain) 
and can be used to start off the parameter continuation. 
Alternatively, the computation can be started from a so- 
lution obtained otherwise (e.g., via time stepping). As 
shown in figure [3j our method automatically generates 
the energy landscape, parametrized by /i, revealing the 
existence of several branches with different energy (for a 
definition of Gibbs free energy see, e.g., [301 Hlj). Ini- 
tially, at low n, the superconducting order parameter is 
strongly suppressed only around the embedded dipole. 
For increasing /i, vortex loops emerge from this area, 
connecting the poles of the dipole. Initially, exactly 4 
such loops are present in the stable solution which en- 
joys the fourfold symmetry of the problem (similar to 
figure 4d). A bifurcation occurs at fj, ps 10.2, suppressing 



three loops towards the dipole and generating a branch of 
states with just one loop (see supplementary material for 
the animation, figure 4a and branch a in figure [3]). Simi- 



lar behavior is found for the state with two loops, with a 
bifurcation point at fi 10.9 (see the state in figure 4b 
and its corresponding energy branch). However, none of 
these states reaches the ground state of the system: it is 
actually the three-loops state that prevails at /i s» 11.4, 
formed via shrinking one loop along the saddle point and 
growing three loops until they hit the cube sides. When 
H is further increased, these loops follow another saddle 
point, where they hit the top and bottom surfaces, and 
then stabilize as three vortices piercing the sample top to 
bottom (figure 4c I. As we continue the computation for 



higher values of /i, more bifurcations occur on each energy 
branch, and new solutions branches emerge. Since this 
calculation is only intended to illustrate our approach, we 
do not include details about such branches. The entire 
solution curves can be obtained from [50] , 

Conclusion and outlook In this letter, we developed 
a computational tool that allows the efficient exploration 
of the steady-state landscape of nonlinear Schrddinger- 
type equations on a high-resolution three-dimensional 
grid. It uses a preconditioned Newton-Krylov solver 
in combination with a numerical parameter continua- 
tion method. The main advantage is that the compu- 
tational cost increases only linearly with the number of 
grid points in the calculation. This is due to the fact 
that the number of required iterations is independent of 
the dimension of the solution space, i.e., the number of 
unknowns. As a result, it is now possible to efficiently 
study three-dimensional physical systems described by 
nonlinear Schrddinger equations even on low-end work- 
stations. Note that our solver is entirely built of existing 
open-source components. In particular we used the high- 
performance continuation solver implemented in LOCA 
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FIG. 3: Numerical continuation in the parameter [i for 
the cube with cavity as depicted in figure [l] The Gibbs 
energy diagram elucidates four solution branches 

spawning from one base branch. The branch \ \ 

corresponds to the superconducting Meissner state and 

contains (ipo;0) when continued to smaller \x. 
Representative solutions on each branch are shown in 
figure |4j 





(c) 



(d) 



FIG. 4: Exemplary states from figure [3] Shown are the 

isosurfaces of |?/>| 2 = 0.1 and the phase of the 
superconducting order parameter on three of the faces 
of the domain. The field |V>| 2 is strongly depleted near 
the magnetic dipole, at the center, where vortex loops 
emerge and stretch out to the boundary. For details, see 
animated data in the supplementary material. 
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By using Newton's method, our approach gives insight 
into saddle point states which are essential for under- 
standing the dynamics of the GL system. At the same 
time, the finite- volume discretization can be used on sam- 
ples with arbitrary shape. We have employed our solver 
for a challenging numerical problem, computing several 
stable and saddle point vortex(-loop) configurations in a 
three-dimensional superconducting cube encapsulating a 
magnetic dipole |32) . 

In conclusion, our method gives access to the dynamics 
of systems that were to date perceived to be too complex 
to be tackled with numerical continuation. It is inher- 
ently applicable to all physical systems modeled by non- 
linear Schrodinger equations ([l}, including the Gross- 
Pitaevskii equations for Bose-Einstein condensates, non- 
linear optics, plasma physics, deep water waves, and even 
seemingly distant subjects such as cosmology and parti- 
cle physics. Just as demonstrated for vortices in super- 
conductors, numerical continuation methods can be used 
for systematic studies of other topological solitons (even 
three-dimensional knotted ones [33J 131]), solitary waves 
and breathers, which are common to various equations 
of nonlinear Schrodinger type (see, e.g., (35|). 
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